%%
%This is the paper FIGURE1, BY YUAN 2021/4/25/ with passion!@#$%
%the section of density and salinity/ along-isopycnal or along-isobath
%curret 
clc;clear;
path(path,'/Users/yuannan/Desktop/sea ice in the okhotsk/in japan/gsw_matlab_v3_06_11');
name_file={'with-amur-5-ave'};
name_month={'Jan.'}; 

eval(['cd /Volumes/Samsung_T5/results_',name_file{1},'/data']);
r1=[130:300];r2=[305:360];
name_var={'s','t','u','v'};
name_var1={'s1','t1','u1','v1'};
name_var2={'s3','t3','u3','v3'};

%read data
for i=1:4
    eval(['fid',num2str(i),'=fopen(','''',name_var{i},'.grd','''',',','''','r','''',');']);
    %eval(['fseek(fid',num2str(i),',800*800*84*4*12*9',',-1);']); %for no fresh 查第10年数据
    eval([name_var{i},'=fread(fid',num2str(i),',','800*800*84',',','''','real*4','''',');']);
    eval([name_var1{i},'=reshape(',name_var{i},',800,800,84);']);
    eval([name_var1{i},'(',name_var1{i},'<-10000)=nan;']);
    eval(['fclose(fid',num2str(i),');']);
    eval([name_var2{i},'=',name_var1{i},'(r1,r2,:);']);
end

%transfer to 56*171*84
p(1:171,1:56,1:84)=10.1325;
rho=gsw_rho(s3,t3,p)-1000;
r_no=permute(rho,[2,1,3]);
s_no=permute(s3,[2,1,3]);
u_no=permute(u3,[2,1,3]);
v_no=permute(v3,[2,1,3]);

% %page inverse
% for z=1:84
%     r_no_1(:,:,z)=r_no(:,:,85-z);
%     s_no_1(:,:,z)=s_no(:,:,85-z);
%     u_no_1(:,:,z)=u_no(:,:,85-z);
%     v_no_1(:,:,z)=v_no(:,:,85-z);
% end

for z=1:84
    r_no_1(:,:,z)=r_no(:,:,z);
    s_no_1(:,:,z)=s_no(:,:,z);
    u_no_1(:,:,z)=u_no(:,:,z);
    v_no_1(:,:,z)=v_no(:,:,z);
end
cd /Volumes/Samsung_T5/program
read_parameter
%cross/along isobath current
for z=1:84
   along(:,:,z)=u_no_1(:,:,z).*ax_new(:,:)+v_no_1(:,:,z).*ay_new(:,:);
end

%coastal thermal wind
dens_sum_1(1:56,1:171)=0;
dens_sum_2(1:56,1:171)=0;
num_dens(1:56,1:171)=1;
dens_ave_1(1:56,1:171)=0;

for i=1:56
    for j=1:171
            for kk=2:84
                if isnan(r_no_1(i,j,kk))
                    continue
                else
                    dens_sum_2(i,j)=(r_no_1(i,j,kk)+r_no_1(i,j,kk-1))./2.0.*(depth(1,kk)-depth(1,kk-1));
                    dens_sum_1(i,j)=dens_sum_1(i,j)+dens_sum_2(i,j);
                    num_dens(i,j)=num_dens(i,j)+1;
                end
            end
        dens_ave_1(i,j)=dens_sum_1(i,j)./depth(1,num_dens(i,j));
    end
end

[densgx,densgy]=gradient(dens_ave_1);
densgx=densgx./x_dis;densgy=densgy./y_dis;
length_n=sqrt(densgx.*densgx+densgy.*densgy);
densgx=densgx./length_n;densgy=densgy./length_n; %gradient direction small--->large; right%downward:positive gradient
densay=sqrt(1./((densgy.*densgy)./(densgx.*densgx)+1));%tangent direction
densax=(-1.*densay.*densgy)./densgx;


%section 53.5 value
for j=1:num_size1
    r_no_lat(j,:)=r_no_1(test3(j,1),test4(j,1),:);
    s_no_lat(j,:)=s_no_1(test3(j,1),test4(j,1),:);
    u_no_lat(j,:)=u_no_1(test3(j,1),test4(j,1),:);
    v_no_lat(j,:)=v_no_1(test3(j,1),test4(j,1),:);
    along_lat(j,:)=along(test3(j,1),test4(j,1),:);
    ax_lat(j,1)=ax_new(test3(j,1),test4(j,1));  
    ay_lat(j,1)=ay_new(test3(j,1),test4(j,1)); 
    densax_lat(j,1)=densax(test3(j,1),test4(j,1));  
    densay_lat(j,1)=densay(test3(j,1),test4(j,1));  
end
%scale_1=[34:41];%53N
scale_1=[39:43];%53.5N
%scale_1=[1:num_size1];
ax_lat(scale_1,1)=densax_lat(scale_1,1);
ay_lat(scale_1,1)=densay_lat(scale_1,1);

for z=1:84
    along_new(:,z)=u_no_lat(:,z).*ax_lat(:,1)+v_no_lat(:,z).*ay_lat(:,1);
end
%%
cd /Volumes/Samsung_T5/paper_figures
fig_scale=[39:52];
depth_c=25;%13:90m 24:200m 19:154m
[off_intep,depth_intep]=meshgrid(off1_dis(2:15,1),depth(1,1:depth_c));
depth_intep(1,:)=0;
fig_num1=2;fig_num2=1;
color_1=[217 39 130]/255;

% h_sub(1)=subplot(fig_num1,fig_num2,1);
% pcolor(off_intep,depth_intep,s_no_lat(fig_scale,1:depth_c)');shading interp;colormap(h_sub(1),flipud(cmocean('haline')));
% hold on
% [cc3,cch3]=contour(off_intep,depth_intep,r_no_lat(fig_scale,1:depth_c)',[22.2 22.4 22.6 22.8 23 23.2 23.4 23.6 23.8 24 24.2 24.4 24.6 24.8 25 25.2 25.4 25.6 25.8 26.0 26.2],...
%     'linewidth',2.5,'color',color_1,'showtext','off');
% hold on;
% [cc3,cch3]=contour(off_intep,depth_intep,r_no_lat(fig_scale,1:depth_c)',[25.7 25.9 26.1],...
%     'linewidth',2.5,'color',color_1,'linestyle',':');
% col3=colorbar;
% caxis([29 33]);set(col3,'ticks',[29 30 31 32 33],'ticklabels',{'29','30','31','32','33'});
% set(col3,'location','southoutside','position',[0.2 0.64 0.35 0.04])
% title(col3,'salinity(psu)','fontsize',20,'FontName','Arial Rounded MT Bold','color','k');%cbfreeze(col3);
% clabel(cc3,cch3,'fontsize',10,'labelspacing',200,'FontName','Arial Black','color',color_1);
% %hold on;text(18,120,'\bf{\rho}','fontsize',30,'FontName','Arial Rounded MT Bold','color','k');
% view(0,-90);

h_sub(2)=figure(1);
pcolor(off_intep,depth_intep,along_new(fig_scale,1:depth_c)');shading interp;colormap(h_sub(2),cmocean('curl'));
hold on
[cc3,cch3]=contour(off_intep,depth_intep,r_no_lat(fig_scale,1:depth_c)',[22 22.2 22.4 22.6 22.8 23 23.2 23.4 23.6 23.8 24 24.2 24.4 24.6 24.8 25 25.2 25.4 25.6 25.8 26.2],...
    'linewidth',2.5,'color',color_1,'showtext','off');
hold on;
[cc3,cch3]=contour(off_intep,depth_intep,r_no_lat(fig_scale,1:depth_c)',[25.7 25.9 26.1],...
    'linewidth',2.5,'color',color_1,'linestyle',':');
clabel(cc3,cch3,'fontsize',15,'labelspacing',200,'FontName','Arial Black','color',color_1);

hold on;
[cc4,cch4]=contour(off_intep,depth_intep,r_no_lat(fig_scale,1:depth_c)',[26 26],...
    'linewidth',2.5,'color',color_1);
clabel(cc4,cch4,'fontsize',15,'labelspacing',200,'FontName','Arial Black','color',color_1);

hold on;
[cc2,cch2]=contour(off_intep,depth_intep,along_new(fig_scale,1:depth_c)',[-25 -20 -15 -10 -5 -1 0],'showtext','on','linewidth',2,'color','k');
clabel(cc2,cch2,'fontsize',15,'labelspacing',200,'FontName','Arial Black','color','k');
col2=colorbar;
%caxis([-10 10]);
caxis([-30 30]);
title(col2,'v(cm/s)','fontsize',20,'FontName','Arial Rounded MT Bold','color','k');
set(col2,'ticks',[-30 0 30],'ticklabels',{'-30','0','30'});
set(col2,'location','southoutside','position',[0.200 0.23 0.35 0.1])
%hold on;text(18,120,'u_a','fontsize',30,'FontName','Arial Rounded MT Bold','color','k');
view(0,-90);
%set(h_sub(2),'outerposition',[0.49 0 0.4281 1]);
set(gca,...
    'ytick',[0 50 100 150 200],...
    'yticklabel',{'0','50','100','150','200'},...
    'xtick',[0 10 20 30 40 50 60],...
    'xticklabel',{'0','10','20','30','40','50','60'},...
    'FontName','Arial Black','fontsize',15,'linewidth',5,...
    'Tickdir','out');
    

set(gcf,'unit','centimeters','position',[0,0,21,13]);
print('ctrl_section1_new','-dtiff','-r600');
